
      MODULE SURFPROP

!SURFACE CHARACTERISITICS FOR DIFFERENT LAND USE CATHEGORIES
!0 - 
!1 - water 
!2 -
!3 -
!4 -
!5 - 
!6 -
      real(4), dimension (0:6)::veg_,wr_,iclay_,w2_,wg_,xirr_,xlai_,top_ &
      &   ,alb_,isa_,z0_,ive_,xdd_,ts_,t2_,rsm_,tsw_,z0w_,tsm_

      contains
      SUBROUTINE SURFPROP_INIT

!0: STANDARD SURFACE

      iclay_(0)=30.
      veg_(0)=0.5
      alb_(0)=0.2
      z0_(0)=0.2
      z0w_(0)=0.01
      ts_(0)=263.
      tsm_(0)=263.
      tsw_(0)=263.
      t2_(0)=263.
      wr_(0)=0.
      wg_(0)=.25
      w2_(0)=.25
      rsm_(0)=30
      xirr_(0)=1.
      ive_(0)=1
      xdd_(0)=.8

!1: BARE SOIL

      iclay_(1)=30.
      veg_(1)=0.
      alb_(1)=0.2
      z0_(1)=0.01
      z0w_(1)=0.01
      ts_(1)=263.
      tsm_(1)=263.
      tsw_(1)=263.
      t2_(1)=263.
      wr_(1)=0.
      wg_(1)=.15
      w2_(1)=.15
      rsm_(1)=30.
      xirr_(1)=1.
      ive_(1)=1
      xdd_(1)=.8

!2: FOREST

      iclay_(2)=30.
      veg_(2)=1.
      alb_(2)=0.2
      z0_(2)=0.2
      z0w_(2)=0.01
      ts_(2)=263.
      tsm_(2)=263.
      tsw_(2)=263.
      t2_(2)=263.
      wr_(2)=0.
      wg_(2)=.30
      w2_(2)=.30
      rsm_(2)=30
      xirr_(2)=1.
      ive_(2)=1
      xdd_(2)=.8

      iclay_(3)=30.
      veg_(3)=0.5
      alb_(3)=0.2
      z0_(3)=0.2
      z0w_(3)=0.01
      ts_(3)=263.
      tsm_(3)=263.
      tsw_(3)=263.
      t2_(3)=263.
      wr_(3)=0.
      wg_(3)=.15
      w2_(3)=.15
      rsm_(3)=30
      xirr_(3)=1.
      ive_(3)=1
      xdd_(3)=.8

      iclay_(4)=30.
      veg_(4)=0.5
      alb_(4)=0.2
      z0_(4)=0.2
      z0w_(4)=0.01
      ts_(4)=263.
      tsm_(4)=263.
      tsw_(4)=263.
      t2_(4)=263.
      wr_(4)=0.
      wg_(4)=.15
      w2_(4)=.15
      rsm_(4)=30
      xirr_(4)=1.
      ive_(4)=1
      xdd_(4)=.8

      iclay_(5)=30.
      veg_(5)=0.5
      alb_(5)=0.2
      z0_(5)=0.2
      z0w_(5)=0.01
      ts_(5)=263.
      tsm_(5)=263.
      tsw_(5)=263.
      t2_(5)=263.
      wr_(5)=0.
      wg_(5)=.15
      w2_(5)=.15
      rsm_(5)=30
      xirr_(5)=1.
      ive_(5)=1
      xdd_(5)=.8

      iclay_(6)=30.
      veg_(6)=0.5
      alb_(6)=0.2
      z0_(6)=0.2
      z0w_(6)=0.01
      ts_(6)=263.
      tsm_(6)=263.
      tsw_(6)=263.
      t2_(6)=263.
      wr_(6)=0.
      wg_(6)=.15
      w2_(6)=.15
      rsm_(6)=30
      xirr_(6)=1.
      ive_(6)=1
      xdd_(6)=.8

      END SUBROUTINE surfprop_init
      END MODULE SURFPROP

      PROGRAM MASKSOIL
      
      use SURFPROP
!
! prepares surface grds (including top)
!
      implicit double precision(a-h,o-z)
!      parameter(nx=193,ny=193)
      parameter (nx = 121, ny = 121)
      dimension xlake(0:nx+1,0:ny+1),veg(0:nx+1,0:ny+1) &
      &  ,wr(0:nx+1,0:ny+1),iclay(0:nx+1,0:ny+1) &
      &  ,w2(0:nx+1,0:ny+1),wg(0:nx+1,0:ny+1),xirr(0:nx+1,0:ny+1) &
      &  ,xlai(0:nx+1,0:ny+1),top(0:nx+1,0:ny+1) &
      &  ,alb(0:nx+1,0:ny+1),isa(0:nx+1,0:ny+1) &
      &  ,z0(0:nx+1,0:ny+1),ive(0:nx+1,0:ny+1),xdd(0:nx+1,0:ny+1) &
      &  ,ts(0:nx+1,0:ny+1),t2(0:nx+1,0:ny+1),rsm(0:nx+1,0:ny+1) &
      &  ,tsw(0:nx+1,0:ny+1),z0w(0:nx+1,0:ny+1) &
      &  ,tsm(0:nx+1,0:ny+1),dep(0:nx+1,0:ny+1) &
      &  ,xice(0:nx+1,0:ny+1),snow(0:nx+1,0:ny+1)
      character*30 fnmap,fext,fnunm,fname
      parameter(wgsat=0.33,w2sat=0.35)
      real(4) dx,dy
      real(4) xx,yy,zz
      real(4), dimension(0:6,0:nx+1,0:ny+1):: landuse
      real(4)  landuse_w(0:6,0:2*nx,0:2*ny), top_w(0:2*nx,0:2*ny)
      integer ii,jj,ilake,itop,i0,j0,ellx,elly,surftyp
!TERRAIN_DOMAIN-4KM FILE IS USED:IX=106, IY=124
!	open (111, file='maskHMAO/topography.dat')
      open (122, file='mask/Surgut_water.02m')
!     open (1,file='mask/Surgut_water.dat',status='old')

      fnmap = 'map/grid_map'

      dx = 1.9E+3
      dy = 1.9E+3

      ilake = 1
! ilake=1 earth surface with lakes
! ilake=0 no lakes
      itop = 0
! itop=1 relief is on
! itop=0 no relief (surface height = 0)

!__________________________________
!  surftyp   | surface landscape  |
!____________|____________________|
!     0      |     "standard"     |
!     1      |      bare soil     |
!     2      |       forest       |
!____________|____________________|

      surftyp = 0
      xlake = 0.
      dep = -1.
      top = 0.

      call SURFPROP_INIT
      landuse(:,:,:)=0.
      
      landuse(surftyp,:,:)=1

!READING LANDUSE AND TOPOGRAPHY DISTRIBUTION FROM FILE
      !do while (.not.eof(1))
      ! read (1,*) i,j,top_w(i,j),(landuse_w(ii,i,j),ii=1,6)
      !enddo
      !do i=0,nx+1
      ! do j=0,ny+1
      !  top(i,j)=top_w(i,j)
      !  landuse(:,i,j)=landuse_w(:,i,j)
      ! enddo 
      !enddo
      !xlake(:,:)=landuse(1,:,:)

!ELLIPTICAL LAKE IN THE CENTER
!      if (ilake == 1) then
!        i0=int((nx+1)/2.)
!        j0=int((ny+1)/2.)
!        ellx=int((nx+1)/5.)
!        elly=int((ny+1)/8.)
!        do i=0,nx+1
!          do j=0,ny+1
!            if ((i-i0)**2/float(ellx**2)+(j-j0)**2/float(elly**2)<=1.) then
!               xlake(i,j)=1.
!               dep(i,j)=2.
!            endif
!          enddo
!        enddo
!      endif   
      
!     xlake(18,18)=1.;xlake(18,19)=1.;xlake(19,18)=1.
!	xlake(18,17)=1.;xlake(17,18)=1.
!	dep(18,18)=2.;dep(18,19)=2.;dep(19,18)=2.
!	dep(18,17)=2.;dep(17,18)=2.

!FOUR ELLIPTICAL LAKES IN THE CORNERS
!      ellx=int((nx+1)/4.)
!	elly=int((ny+1)/4.)
!      do i=0,nx+1
!	 do j=0,ny+1
!        if ( (i-1)**2/float(ellx**2)+(j-1)**2/float(elly**2)<=1.or.
!     &      (i-1)**2/float(ellx**2)+(j-ny)**2/float(elly**2)<=1.or.
!     &	  (i-nx)**2/float(ellx**2)+(j-1)**2/float(elly**2)<=1.or.
!     &	  (i-nx)**2/float(ellx**2)+(j-ny)**2/float(elly**2)<=1.
!     &    ) then
!	   xlake(i,j)=1.
!	   dep(i,j)=2.
!        endif	  	
!	 enddo
!      enddo

!	xlake=0.

!T-FORM WATER BODY DISTRIBUTION (A RIVER AND A TRIBUTARY)
!ir=nx/2
!jr=ny/2
!do j = 0,ny+1
!   xlake(ir,j)=0.
      !end do
!do i=ir+1,nx+1
!   xlake(i,jr)=0.
!end do

!SEA (EASTERN PART) AND LAND (WESTERN PART)
      !if (ilake==1) then
! do i = 1, int(float(nx)/2.)
!  do j = 1, ny
!   xlake(i,j)=1.
      !  enddo
      ! enddo
      !endif

!16 RECTANGULAR LAKES
     
!      if (ilake==1) then
!	 do ii=1,7,2
!	  do i=ii*4,(ii+1)*4
!	   do jj = 1,7,2
!	    do j=jj*4,(jj+1)*4
!	     xlake(i,j)=1
!	     dep(i,j)=10.
!          enddo
!	   enddo
!       enddo
!	 enddo   
!      endif

!RECTANGLE SHAPED LAKE IN CENTER OF DOMAIN
      !if (ilake==1) then
      ! do i=int(float(nx)/4.), int(3.*float(nx)/4.)
!  do j=int(float(ny)/4.), int(3.*float(ny)/4.)
!   xlake(i,j)=1.
      !  enddo
      ! enddo
!endif

!READING FILE WITH WATER BODY DISTRIBUTION 
!xlake(i,j): 0 - soil, 1 - water
!      if (ilake==1) then
!       do j = 1, 60 !102
!	  do i = 1, 120 !123
!	   if ((i>=1.and.i<=2*(nx+1)+1).and.(j>=1.and.j<=ny+2).and.
!     &   	 mod(i-1,2)==0) then
!	    read(122,*) ii, jj, xlake((i-1)/2,ny+2-j)
!	    dep((i-1)/2,ny+2-j) = 2.
!         else
!	    read(122,*) ii, jj, xx
!         endif
!        enddo
!       enddo
!      endif 
!      xlake = xlake/100.	

!      if (ilake == 1) then
!        do j = 1, 60 !102
!          do i = 1, 120 !123
!            read(122,*) xx, yy, zz
!            if (j == 1) then
!              xlake(i-1,119) = zz
!              dep(i-1,119) = 2.
!            else  
!              xlake(i-1,120-2*j+2) = 0.5*(zz + xlake(i-1,120-2*j+2))
!              xlake(i-1,120-2*j+1) = zz
!              dep(i-1,120-2*j+1:120-2*j+2) = 2.
!            endif
!          enddo
!        enddo
!      endif
      
!      do i = 0, 119
!        xlake(i,120:ny+1) = xlake(i,119)
!      enddo
!      do j = 0, ny+1
!        xlake(120:nx+1,j) = xlake(119,j)
!      enddo

!      xlake = xlake/100.


!READING FILE WITH TOPOGRAPHY
!      if (itop==1) then
!       do i=1, 106  !102
!	  do j=1, 124  !123
!	   if ((i>=1.and.i<=1+nx+1).and.(j>=1.and.j<=1+ny+1)) then
!	    read(111,*) ii, jj, top(i-1,j-1)
!         else
!	    read(111,*) ii, jj, xx
!         endif
!        enddo
!       enddo
!      endif
!do i=0,102
! do j=124,127
!  top(i,j) = top(i,123)
      ! enddo
      !enddo
      

      xice(:,:) = 0. ! lake/sea ice depth, m
      snow(:,:) = 0. ! snow depth over land and lake/sea ice, m
      xdd(:,:) = 10. ! "Soil" depth, m
      ts(:,:) = 293. ! Soil temperature, K
      wg(:,:) = 0.35 ! Soil wetness, kg/kg


      call wrigrd(xlake,xmi,xma,ymi,yma,vmi,vma,dx,dy &
      &   ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'lak'))
      call wrigrd(dep,xmi,xma,ymi,yma,vmi,vma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'dep'))

      albs=0.2
      albw=0.05
      do j=0,ny+1
         do i=0,nx+1
            xlai(i,j)=2.
            isa(i,j)=35
!           top=0.
            if (xlake(i,j).eq.1) then
               iclay(i,j)=-2
               veg(i,j)=0.
               alb(i,j)=0.05
               z0w(i,j)=0.01
               z0(i,j)=0.01
               ts(i,j)=273.
               tsm(i,j)=273.
               tsw(i,j)=273.
               t2(i,j)=273.
               wr(i,j)=0.
               wg(i,j)=.15
               w2(i,j)=.15
               rsm(i,j)=30
               xirr(i,j)=0.
               ive(i,j)=1
               xdd(i,j)=.8
             else
               iclay(i,j)=sum(landuse(:,i,j)*iclay_(:))
               veg(i,j)=  sum(landuse(:,i,j)*veg_(:))
               alb(i,j)=  sum(landuse(:,i,j)*alb_(:))
               z0(i,j)=   sum(landuse(:,i,j)*z0_(:))
               z0w(i,j)=  sum(landuse(:,i,j)*z0w_(:))
!               ts(i,j)=   sum(landuse(:,i,j)*ts_(:))
               tsm(i,j)=  sum(landuse(:,i,j)*tsm_(:))
               tsw(i,j)=  sum(landuse(:,i,j)*tsw_(:))
               t2(i,j)=   sum(landuse(:,i,j)*t2_(:))
               wr(i,j)=   sum(landuse(:,i,j)*wr_(:))
!               wg(i,j)=   sum(landuse(:,i,j)*wg_(:))
               w2(i,j)=   sum(landuse(:,i,j)*w2_(:))
               rsm(i,j)=  sum(landuse(:,i,j)*rsm_(:))
               xirr(i,j)=  sum(landuse(:,i,j)*xirr_(:))
               ive(i,j)=  sum(landuse(:,i,j)*ive_(:))
!               xdd(i,j)=  sum(landuse(:,i,j)*xdd_(:))
            endif
         enddo
      enddo

!	16 WET SOIL PATCHES
!	do ii=1,7,2
!	 do i=ii*4,(ii+1)*4
!	  do jj = 1,7,2
!	   do j=jj*4,(jj+1)*4
!	    veg(i,j)=0.
!	    w2(i,j)=.30
!          wg(i,j)=.30
!         enddo
!	  enddo
!       enddo
!	enddo   

      call iwrigrd(iclay,xmi,xma,ymi,yma,icmi,icma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'icl'))
      call wrigrd(veg,xmi,xma,ymi,yma,vmi,vma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'veg'))
      call wrigrd(w2,xmi,xma,ymi,yma,w2mi,w2ma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'w2'))
      call wrigrd(wg,xmi,xma,ymi,yma,wgmi,wgma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'wg'))
      call wrigrd(xlai,xmi,xma,ymi,yma,wgmi,wgma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'lai'))
      call wrigrd(top,xmi,xma,ymi,yma,wgmi,wgma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'top'))
      call wrigrd(alb,xmi,xma,ymi,yma,albmi,albma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'alb'))
      call wrigrd(z0,xmi,xma,ymi,yma,z0mi,z0ma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'z0'))
      call wrigrd(z0w,xmi,xma,ymi,yma,z0wmi,z0wma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'z0w'))
      call wrigrd(ts,xmi,xma,ymi,yma,tsmi,tsma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'ts'))
      call wrigrd(tsm,xmi,xma,ymi,yma,tsmi,tsma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'tsm'))
      call wrigrd(tsw,xmi,xma,ymi,yma,tswmi,tswma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'tsw'))
      call wrigrd(t2,xmi,xma,ymi,yma,t2mi,t2ma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'t2'))
      call wrigrd(wr,xmi,xma,ymi,yma,wrmi,wrma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'wr'))
      call wrigrd(xirr,xmi,xma,ymi,yma,xirrmi,xirrma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'irr'))
      call iwrigrd(isa,xmi,xma,ymi,yma,isami,isama,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'isa'))
      call iwrigrd(ive,xmi,xma,ymi,yma,ivemi,ivema,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'ive'))
      call wrigrd(rsm,xmi,xma,ymi,yma,rsmmi,rsmma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'rsm'))
      call wrigrd(xdd,xmi,xma,ymi,yma,xddmi,xddma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'xdd'))
      call wrigrd(xice,xmi,xma,ymi,yma,xddmi,xddma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'ice'))
       call wrigrd(snow,xmi,xma,ymi,yma,xddmi,xddma,dx,dy &
      &  ,0,nx+1,0,ny+1,0,nx+1,0,ny+1,fext(fnmap,'snw'))

      END PROGRAM MASKSOIL

      subroutine wrigrd(z,xmin,xmax,ymin,ymax,zmin,zmax,dx,dy &
      &   ,m0,m1,n0,n1,i0,i1,j0,j1,fname)

      implicit double precision (a-h,o-z)
      character*30 fname
      dimension z(m0:m1,n0:n1)
      real(4) dx,dy

      zmin=z(i0,j0)
      zmax=z(i0,j0)

      do j=j0,j1
         do i=i0,i1
           zmin=min(zmin,z(i,j))
           zmax=max(zmax,z(i,j))
         enddo
      enddo

      xmin=-float(i1)/2.*dx/1000.
      xmax=float(i1)/2.*dx/1000.
      ymin=-float(j1)/2.*dy/1000.
      ymax=float(j1)/2.*dy/1000.

      open(77,file=fname)
      write(77,'(a5)') 'DSAA '
      write(77,*) (i1-i0+1),(j1-j0+1), ' '
      write(77,*) xmin,xmax, ' '
      write(77,*) ymin,ymax, ' '
      write(77,*) zmin,zmax, ' '

      do j=j0,j1
         write(77,fmt=1) (z(i,j),i=i0,i1)
      enddo

      close(77)
1     format (<i1-i0+1>E15.7, " ")
      return
      end

      subroutine iwrigrd(iz,xmin,xmax,ymin,ymax,izmin,izmax,dx,dy &
      &   ,m0,m1,n0,n1,i0,i1,j0,j1,fname)

      implicit double precision (a-h,o-z)
      character*30 fname
      dimension iz(m0:m1,n0:n1)
      real(4) dx,dy

      izmin=iz(i0,j0)
      izmax=iz(i0,j0)

      do j=j0,j1
         do i=i0,i1
           izmin=min(izmin,iz(i,j))
           izmax=max(izmax,iz(i,j))
         enddo
      enddo

      xmin=-float(i1)/2.*dx/1000.
      xmax=float(i1)/2.*dx/1000.
      ymin=-float(j1)/2.*dy/1000.
      ymax=float(j1)/2.*dy/1000.

      open(77,file=fname)
      write(77,'(a5)') 'DSAA '
      write(77,*) (i1-i0+1),(j1-j0+1), ' '
      write(77,*) xmin,xmax, ' '
      write(77,*) ymin,ymax, ' '
      write(77,*) izmin,izmax, ' '

      do j=j0,j1
         write(77,fmt = 2) (iz(i,j),i=i0,i1)
      enddo

      close(77)
2     format (<i1-i0+1>i8, " ")
      return
      end

      subroutine readgrd(nunit,var,i0,i1,j0,j1,x0,x1,y0,y1,z0,z1)
!
! reads 2d array from a grd file
!
      implicit double precision(a-h,o-z)
      dimension var(i0:i1,j0:j1)
      character*4 dsaa

      read(nunit,'(a4)') dsaa
      read(nunit,*) nx,ny
      read(nunit,*) x0,x1
      read(nunit,*) y0,y1
      read(nunit,*) z0,z1

      do j=j0,j1
         read(nunit,*) (var(i,j),i=i0,i1)
      enddo

      return
      end

      subroutine ireadgrd(nunit,ivar,i0,i1,j0,j1,x0,x1,y0,y1,iz0,iz1)
!
! reads 2d array from a grd file
!
      implicit double precision(a-h,o-z)
      dimension ivar(i0:i1,j0:j1)
      character*4 dsaa

      read(nunit,'(a4)') dsaa
      read(nunit,*) nx,ny
      read(nunit,*) x0,x1
      read(nunit,*) y0,y1
      read(nunit,*) iz0,iz1

      do j=j0,j1
         read(nunit,*) (ivar(i,j),i=i0,i1)
      enddo

      return
      end

      function fext(fname,ext)
      character*30 fname,fext
      character*3 ext
      parameter(n=40)
!
! elimina extensao se existir
!
      fext=fname

      do 10 i=30,1,-1
         if(fext(i:i).eq.'.') then
            fext(i:n)=char(0)
            go to 11
         endif
10    continue
11    continue

      do 20 i=30,1,-1
         if(fext(i:i).ne.' ' .and. fext(i:i).ne.char(0)) then
            fext(i+1:i+4)='.'//ext
            go to 21
         endif
20    continue
      write(*,*) 'erro em fext'
      pause
21    continue
      return
      end
